knitr::opts_chunk$set(echo = TRUE, comment=" ", cache=TRUE, warning=FALSE, message=FALSE, fig.width=14, fig.asp=.678)
library(openxlsx)
library(data.table)
library(RColorBrewer)
pal <- brewer.pal(5,"Set2")
#library(heatmaply)
library(bipartite)
## Loading required package: vegan
## Loading required package: permute
## Loading required package: lattice
## This is vegan 2.6-4
## Loading required package: sna
## Loading required package: statnet.common
##
## Attaching package: 'statnet.common'
## The following objects are masked from 'package:base':
##
## attr, order
## Loading required package: network
##
## 'network' 1.18.2 (2023-12-04), part of the Statnet Project
## * 'news(package="network")' for changes since last version
## * 'citation("network")' for citation information
## * 'https://statnet.org' for help, support, and other information
## sna: Tools for Social Network Analysis
## Version 2.7-2 created on 2023-12-05.
## copyright (c) 2005, Carter T. Butts, University of California-Irvine
## For citation information, type citation("sna").
## Type help(package="sna") to get started.
## This is bipartite 2.19.
## For latest changes see versionlog in ?"bipartite-package". For citation see: citation("bipartite").
## Have a nice time plotting and analysing two-mode networks.
##
## Attaching package: 'bipartite'
## The following object is masked from 'package:vegan':
##
## nullmodel
data(mosquin1967)
library(vegan)
library(phytools)
## Loading required package: ape
##
## Attaching package: 'ape'
## The following objects are masked from 'package:sna':
##
## consensus, degree
## Loading required package: maps
##
## Attaching package: 'phytools'
## The following object is masked from 'package:vegan':
##
## scores
library(taxize)
#remotes::install_github("ropensci/brranching")
library(brranching)
library(rotl)
##
## Attaching package: 'rotl'
## The following objects are masked from 'package:taxize':
##
## synonyms, tax_name, tax_rank
library(dendextend)
## Registered S3 method overwritten by 'dendextend':
## method from
## rev.hclust vegan
##
## ---------------------
## Welcome to dendextend version 1.17.1
## Type citation('dendextend') for how to cite the package.
##
## Type browseVignettes(package = 'dendextend') for the package vignette.
## The github page is: https://github.com/talgalili/dendextend/
##
## Suggestions and bug-reports can be submitted at: https://github.com/talgalili/dendextend/issues
## You may ask questions at stackoverflow, use the r and dendextend tags:
## https://stackoverflow.com/questions/tagged/dendextend
##
## To suppress this message use: suppressPackageStartupMessages(library(dendextend))
## ---------------------
##
## Attaching package: 'dendextend'
## The following object is masked from 'package:phytools':
##
## untangle
## The following objects are masked from 'package:ape':
##
## ladderize, rotate
## The following object is masked from 'package:permute':
##
## shuffle
## The following object is masked from 'package:data.table':
##
## set
## The following object is masked from 'package:stats':
##
## cutree
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.0 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::between() masks data.table::between()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::first() masks data.table::first()
## ✖ lubridate::hour() masks data.table::hour()
## ✖ lubridate::isoweek() masks data.table::isoweek()
## ✖ dplyr::lag() masks stats::lag()
## ✖ dplyr::last() masks data.table::last()
## ✖ purrr::map() masks maps::map()
## ✖ lubridate::mday() masks data.table::mday()
## ✖ lubridate::minute() masks data.table::minute()
## ✖ lubridate::month() masks data.table::month()
## ✖ lubridate::quarter() masks data.table::quarter()
## ✖ lubridate::second() masks data.table::second()
## ✖ purrr::transpose() masks data.table::transpose()
## ✖ lubridate::wday() masks data.table::wday()
## ✖ lubridate::week() masks data.table::week()
## ✖ dplyr::where() masks ape::where()
## ✖ lubridate::yday() masks data.table::yday()
## ✖ lubridate::year() masks data.table::year()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(heatmaply)
## Loading required package: plotly
##
## Attaching package: 'plotly'
##
## The following object is masked from 'package:ggplot2':
##
## last_plot
##
## The following object is masked from 'package:stats':
##
## filter
##
## The following object is masked from 'package:graphics':
##
## layout
##
## Loading required package: viridis
## Loading required package: viridisLite
##
## Attaching package: 'viridis'
##
## The following object is masked from 'package:maps':
##
## unemp
##
## Registered S3 method overwritten by 'seriation':
## method from
## reorder.hclust vegan
##
## ======================
## Welcome to heatmaply version 1.5.0
##
## Type citation('heatmaply') for how to cite the package.
## Type ?heatmaply for the main documentation.
##
## The github page is: https://github.com/talgalili/heatmaply/
## Please submit your suggestions and bug-reports at: https://github.com/talgalili/heatmaply/issues
## You may ask questions at stackoverflow, use the r and heatmaply tags:
## https://stackoverflow.com/questions/tagged/heatmaply
## ======================
library(googlesheets4)
#load network
#load network
setwd("~/Desktop/Scent-Net")
gs4_auth(email=TRUE)
net.long <- read_sheet("1plhxbVh5IQLNVMazqO4OvfM3X6vbll_Y6BTMJ94bkI4", na=c("","na", "NA")) %>%
mutate(pollinator_group= if_else(pollinator=="papilio_gothica", "butterfly", pollinator_group)) %>% #fix papilio_gothica from bee -> butterfly
mutate(plant= recode(plant,senecio_tall_int = "senecio_integerrimus", agoceris_glauca = "agoseris_glauca", agoceris_aurantiaca="agoseris_aurantiaca", hydrophyllum_fendeleri="hydrophyllum_fendleri"))
net <- net.long %>% select(plant, pollinator, interactions) %>% pivot_wider(names_from = plant, values_from = interactions, values_fn= sum, values_fill= 0) %>% #pivots data set longer so plants as columns and pollinators as rows, replaces NAs with zeros
column_to_rownames("pollinator") #pollinators are in first column, makes them into rownames; #deletes the first pollinators column
net <- as.matrix(t(net)) # turns it from dataframe to matrix bc packages like matrix; t transposes (turns dataset sideways - rows -> columns; columns -> rows)
polls <- colnames(net) #now pollinators are column names
plants<-rownames(net)
snet <- sortweb(net) #sort by row & column totals; sorted rows by rowsums and columns by column sums. Biggest numbers in top left of matrix
#load scent
#Go to comm_vols for code to make scent.net
#Loads in vol, filters with bouquet, takes mean of each voc for each species
#From comm_vols.Rmd
scent.net <- read_csv("data/Volatiles_by_species_mean.csv") %>% column_to_rownames("species") %>% #species are now rownames
as.matrix() #network analysis likes matrix
scent.net <- scent.net[,colnames(scent.net)!="Cyclohexane, 1,3,5-trimethyl-2-octadecyl-"] #filtering this compound bc it has 0 interactions or emissions... 355 vs 354. breaks code later
#loads in chemical class and gets rid of ?
compound.class <-read_sheet("1BlDIOs4STe5ZTTPYE6qez1Nsey_ovEpe5m-SN8Q55gg", sheet= "chemical_class") %>% mutate(major_class= factor(major_class,levels=c("aliphatic", "monoterpene", "sesquiterpene", "nitrogenous", "sulfur", "benzenoid")) %>% fct_na_value_to_level("other")) %>%
filter(compound_long!="Cyclohexane, 1,3,5-trimethyl-2-octadecyl-") #filtering this compound bc it has 0 interactions or emissions... 355 vs 354
#codes chemical class by colors
voc.class.colors <- setNames(brewer.pal(7,"Set3"), levels(compound.class$major_class))
##match paul’s plants with the plants we have vocs for
scentplants<-rownames(scent.net) #%>% tolower() %>% str_replace(" ", "_")
paulsplants <- plants
#plants not in both datasets
setdiff(scentplants, paulsplants) #gives us bonus plants
[1] "aquilegia_caerulea" "frasera_speciosa" "geum_aleppicum"
setdiff(paulsplants, scentplants)
[1] "hydrophyllum_fendleri" "phacelia_heterophylla"
bothplants <- intersect(scentplants, paulsplants) #plants found in amanda and paul data
polls.count <-count(net.long, pollinator_group,pollinator_family, pollinator) #counting number of rows for each pollinator
polls.class <- tibble(pollinator=polls) %>% left_join(polls.count) %>% pull(pollinator_group)
polls.class.dataframe <- tibble(pollinator=polls) %>% left_join(polls.count)
polls.family <- tibble(pollinator=polls) %>% left_join(polls.count) %>% pull(pollinator_family)
polls.order.all <- polls.class
#net.grp.all has all of pauls 45 plant species
#net.grp only has plant species found in Amanda's and paul's datasets
net.grp.all <- aggregate(.~ polls.class, as.data.frame(t(net)), sum) #lump network by pollinator class; plants columns, pollinators rows
rownames(net.grp.all) <- net.grp.all[,1] #make plants rows, pollinators column names
net.grp.all[,1] <- NULL
net.grp.all <- as.matrix(t(net.grp.all)) #transposes matrix
net.grp<- net.grp.all[bothplants,]
net.cut <- net[bothplants,]
##Pollinator groups
#Tells you plants major pollinator
plants.majorpoll <- setNames(factor(colnames(net.grp)[apply(net.grp, 1, which.max)]), rownames(net.grp))
#sets colors for ALL 5 pollinator groups
pcols.all <- setNames(brewer.pal(5,"Set2"), colnames(net.grp))
#Four colors for each major pollinator group
pcols <- pcols.all[ levels(plants.majorpoll)]
#gives % bee pollination for each plant species
plants.bees <- net.grp[,"bee"]/rowSums(net.grp)
plot(sort(plants.bees))

scent.net.cut <- scent.net[rownames(scent.net)%in% paulsplants, ] #cuts scent to match what paul has pollinator data for
# Cap of plants.bees
plant.bee.cap <- capscale(sqrt(scent.net.cut)~plants.bees)
anova(plant.bee.cap)
Permutation test for capscale under reduced model
Permutation: free
Number of permutations: 999
Model: capscale(formula = sqrt(scent.net.cut) ~ plants.bees)
Df Variance F Pr(>F)
Model 1 505636 1.0525 0.35
Residual 41 19696248
#p=0.348
plot(plant.bee.cap, type="text")

#gives % fly pollination for each plant species
plants.fly <- net.grp[,"fly"]/rowSums(net.grp)
plot(sort(plants.fly))

sort(plants.fly)
agoseris_aurantiaca amelanchier_alnifolia
0.000000000 0.000000000
castilleja_sulphurea gentiana_parryi
0.000000000 0.000000000
ipomopsis_aggregata lathyrus_lanszwertii
0.000000000 0.000000000
mahonia_repens mertensia_ciliata
0.000000000 0.000000000
ranunculus_inamoenus rosa_woodsii
0.000000000 0.000000000
vicia_americana viola_praemorsa
0.000000000 0.000000000
delphinium_nuttallianum hydrophyllum_capitatum
0.001774623 0.010204082
campanula_rotundifolia mertensia_fusiformis
0.011428571 0.011627907
erythronium_grandiflorum heliomeris_multiflora
0.019230769 0.033676658
sedum_lanceolatum heterotheca_villosa
0.035398230 0.044919786
lomatium_dissectum taraxacum_officinale
0.076923077 0.102841678
senecio_serra senecio_integerrimus
0.125000000 0.174044876
arenaria_congesta boechera_stricta
0.178947368 0.181818182
erigeron_speciosus solidago_multiradiata
0.202226345 0.205357143
claytonia_lanceolata helianthella_quinquenervis
0.244186047 0.267904509
fragaria_virginiana linum_lewisii
0.333333333 0.333333333
potentilla_hippiana erigeron_flagellaris
0.397425583 0.432900433
potentilla_gracilis eriogonum_umbellatum
0.433962264 0.442477876
pseudocymopterus_montanus draba_aurea
0.490196078 0.564285714
eriogonum_subalpinum androsace_septentrionalis
0.575000000 0.630000000
galium_boreale achillea_millefolium
0.875000000 0.942622951
agoseris_glauca
1.000000000
#gives % butterfly pollination for each plant species
plants.butterfly <- net.grp[,"butterfly"]/rowSums(net.grp)
plot(sort(plants.butterfly))

sort(plants.butterfly)
agoseris_glauca amelanchier_alnifolia
0.000000000 0.000000000
androsace_septentrionalis campanula_rotundifolia
0.000000000 0.000000000
castilleja_sulphurea draba_aurea
0.000000000 0.000000000
erythronium_grandiflorum fragaria_virginiana
0.000000000 0.000000000
gentiana_parryi hydrophyllum_capitatum
0.000000000 0.000000000
linum_lewisii lomatium_dissectum
0.000000000 0.000000000
mahonia_repens mertensia_ciliata
0.000000000 0.000000000
mertensia_fusiformis pseudocymopterus_montanus
0.000000000 0.000000000
ranunculus_inamoenus rosa_woodsii
0.000000000 0.000000000
sedum_lanceolatum senecio_serra
0.000000000 0.000000000
vicia_americana viola_praemorsa
0.000000000 0.000000000
heliomeris_multiflora potentilla_hippiana
0.003934730 0.006436042
delphinium_nuttallianum taraxacum_officinale
0.006654836 0.008119080
achillea_millefolium helianthella_quinquenervis
0.012295082 0.013262599
solidago_multiradiata potentilla_gracilis
0.013392857 0.018867925
heterotheca_villosa claytonia_lanceolata
0.022887701 0.023255814
lathyrus_lanszwertii arenaria_congesta
0.025000000 0.057142857
senecio_integerrimus erigeron_flagellaris
0.065494239 0.089466089
eriogonum_subalpinum galium_boreale
0.125000000 0.125000000
eriogonum_umbellatum erigeron_speciosus
0.150442478 0.179962894
ipomopsis_aggregata boechera_stricta
0.217391304 0.363636364
agoseris_aurantiaca
1.000000000
#gives % hummingbird pollination for each plant species
plants.hummingbird <- net.grp[,"hummingbird"]/rowSums(net.grp)
plot(sort(plants.hummingbird))

sort(plants.hummingbird)
achillea_millefolium agoseris_aurantiaca
0.000000000 0.000000000
agoseris_glauca amelanchier_alnifolia
0.000000000 0.000000000
androsace_septentrionalis arenaria_congesta
0.000000000 0.000000000
campanula_rotundifolia castilleja_sulphurea
0.000000000 0.000000000
claytonia_lanceolata draba_aurea
0.000000000 0.000000000
erigeron_flagellaris erigeron_speciosus
0.000000000 0.000000000
eriogonum_subalpinum eriogonum_umbellatum
0.000000000 0.000000000
erythronium_grandiflorum fragaria_virginiana
0.000000000 0.000000000
galium_boreale gentiana_parryi
0.000000000 0.000000000
helianthella_quinquenervis heliomeris_multiflora
0.000000000 0.000000000
heterotheca_villosa hydrophyllum_capitatum
0.000000000 0.000000000
lathyrus_lanszwertii linum_lewisii
0.000000000 0.000000000
lomatium_dissectum mahonia_repens
0.000000000 0.000000000
mertensia_ciliata mertensia_fusiformis
0.000000000 0.000000000
potentilla_gracilis potentilla_hippiana
0.000000000 0.000000000
pseudocymopterus_montanus ranunculus_inamoenus
0.000000000 0.000000000
rosa_woodsii sedum_lanceolatum
0.000000000 0.000000000
senecio_serra solidago_multiradiata
0.000000000 0.000000000
taraxacum_officinale vicia_americana
0.000000000 0.000000000
viola_praemorsa senecio_integerrimus
0.000000000 0.002425713
boechera_stricta delphinium_nuttallianum
0.090909091 0.224933452
ipomopsis_aggregata
0.608695652
#gives % wasp pollination for each plant species
plants.wasp <- net.grp[,"wasp"]/rowSums(net.grp)
plot(sort(plants.wasp))

sort(plants.wasp)
achillea_millefolium agoseris_aurantiaca
0.0000000000 0.0000000000
agoseris_glauca amelanchier_alnifolia
0.0000000000 0.0000000000
androsace_septentrionalis arenaria_congesta
0.0000000000 0.0000000000
boechera_stricta campanula_rotundifolia
0.0000000000 0.0000000000
castilleja_sulphurea claytonia_lanceolata
0.0000000000 0.0000000000
delphinium_nuttallianum draba_aurea
0.0000000000 0.0000000000
erigeron_flagellaris erigeron_speciosus
0.0000000000 0.0000000000
eriogonum_umbellatum erythronium_grandiflorum
0.0000000000 0.0000000000
fragaria_virginiana galium_boreale
0.0000000000 0.0000000000
gentiana_parryi helianthella_quinquenervis
0.0000000000 0.0000000000
heliomeris_multiflora hydrophyllum_capitatum
0.0000000000 0.0000000000
ipomopsis_aggregata lathyrus_lanszwertii
0.0000000000 0.0000000000
linum_lewisii lomatium_dissectum
0.0000000000 0.0000000000
mahonia_repens mertensia_ciliata
0.0000000000 0.0000000000
mertensia_fusiformis potentilla_hippiana
0.0000000000 0.0000000000
pseudocymopterus_montanus ranunculus_inamoenus
0.0000000000 0.0000000000
rosa_woodsii sedum_lanceolatum
0.0000000000 0.0000000000
senecio_integerrimus senecio_serra
0.0000000000 0.0000000000
taraxacum_officinale vicia_americana
0.0000000000 0.0000000000
viola_praemorsa heterotheca_villosa
0.0000000000 0.0005347594
solidago_multiradiata potentilla_gracilis
0.0178571429 0.0377358491
eriogonum_subalpinum
0.0750000000
##Heatmap1
#pp_heatmap
heatmaply(net, scale="none", showticklabels = c(F,T), margins= c(0,NA,0,0),
scale_fill_gradient_fun = ggplot2::scale_fill_gradient(low = "white", high = "#000050", trans="log1p"),
distfun=vegdist,
hclustfun=function(x) hclust(x, method="mcquitty"),
row_side_colors = data.frame(P=plants.majorpoll[rowSums(net)>0]),
row_side_palette=function(n) pcols,
col_side_colors=data.frame(PollinatorGroup=polls.class),
col_side_palette=function(n) pcols.all)
##Matrix plot
#interaction matrix plot
visweb(net, type="nested", circles=T, boxes=F, circle.max=1.2)

##Web plot
#fix names of species (plants rows, polls cols)
net.nice.labels <- net
rownames(net.nice.labels) <- str_to_sentence(str_replace_all(rownames(net), "_"," "))
colnames(net.nice.labels) <- str_to_sentence(str_replace_all(colnames(net), "_"," "))
#interaction web plot
par(mar=c(0,0,0,0))
plant.pollinator.network <- plotweb(sqrt(net.nice.labels), empty=F, text.rot=90, method="cca", col.high=pcols.all[polls.order.all], col.low=pcols[plants.majorpoll], col.interaction=pcols.all[polls.order.all], bor.col.interaction=pcols.all[polls.order.all], bor.col.high=pcols.all[polls.order.all], bor.col.low=pcols[plants.majorpoll])

library(ggplot2)
ggsave(plant.pollinator.network, file='plant.pollinator.networkplot.jpg', width=130, height=105, units= 'mm', dpi=600)
#computes network statistics (like nestedness)
#networkstats <- networklevel(net)
#write_csv(enframe(networkstats), "data/p-p network stats.csv")
#computes network statistics for voc-pollinator network
#scentnetworkstats<- networklevel(scentpollnet)
#write_csv(enframe(scentnetworkstats), "data/scent network stats.csv")
scentnet.stats <- read_csv("data/scent network stats.csv")
#the network statistics null model: TODO- broken
#network.nullmodel<-nullmodel(scentpollnet)
#nullmodel.stats<- map(network.nullmodel, ~networklevel(.x )) #network.nullmodel[1:2]; index=c("NODF", "weighted NODF", "connectance", "H2", "modularity Q")
#nullmodel.stats %>% map_dfr(~enframe(.x), .id = "run") %>% write_csv("data/scent network null model statistics.csv")
nullmodel.stats <- read_csv("data/scent network null model statistics.csv")
#idea - pivot wider so network properties are the columns, rows are run number, then average by values in column -> pivot longer or store in new variable
#nullmodel.stats %>% pivot_wider(names_from = "Name", values_from = "Value") %>% column_to_rownames("stats") %>% as.matrix()
null.NODF <- read_csv("data/scent network null model statistics.csv") %>% filter(name=="NODF")
ggplot(null.NODF, aes(x=value))+ geom_histogram() + geom_vline(aes(xintercept=filter(scentnet.stats, name=="NODF") %>% pull(value)))

##Heatmap2
scent.net.cut <- scent.net[rownames(scent.net)%in% paulsplants, ] #cuts scent to match what paul has pollinator data for
plants.majorpoll.cut <- plants.majorpoll[rownames(scent.net.cut)] #cuts Pauls data to match Amanda's
plants.bees.cut<- plants.bees[rownames(scent.net.cut)]
heatmaply(sqrt(scent.net.cut), scale="none", showticklabels = c(F,T), margins= c(0,NA,0,0),
scale_fill_gradient_fun = ggplot2::scale_fill_gradient(low = "white", high = "#000050", trans="log1p"),
distfun=vegdist,
hclustfun=function(x) hclust(x, method="mcquitty"),
row_side_colors = data.frame(major.polls=plants.majorpoll.cut),
row_side_palette=function(n) pcols )
#col_side_colors=data.frame(PollinatorGroup=scent.net), #put compound class here
# col_side_palette=function(n) pcols.all)
#scent-plant network
distance.matrix <- vegdist(scent.net.cut)
View(as.matrix(distance.matrix)) #tells you how different the smells are. If 0 there is no difference, scents are the same; 1 = no compounds in common
#plant-pollinator network
#how similar are the plants in terms of their pollinator array/composition
distance.matrix.pp <- vegdist(net.cut)
View(as.matrix(distance.matrix.pp))
#Mantel Test
#can do correlation, if scent similarity predicts pollinator --> mantel test
mantel(distance.matrix,distance.matrix.pp)
Mantel statistic based on Pearson's product-moment correlation
Call:
mantel(xdis = distance.matrix, ydis = distance.matrix.pp)
Mantel statistic r: 0.02341
Significance: 0.321
Upper quantiles of permutations (null model):
90% 95% 97.5% 99%
0.0690 0.0835 0.0971 0.1197
Permutation: free
Number of permutations: 999
#TODO: group by family then rerun vegdist
distance.matrix.pollgroup <- vegdist(net.grp)
mantel(distance.matrix, distance.matrix.pollgroup)
Mantel statistic based on Pearson's product-moment correlation
Call:
mantel(xdis = distance.matrix, ydis = distance.matrix.pollgroup)
Mantel statistic r: -0.06377
Significance: 0.892
Upper quantiles of permutations (null model):
90% 95% 97.5% 99%
0.0670 0.0922 0.1132 0.1510
Permutation: free
Number of permutations: 999
#can't predict which pollinator group will visit plant based on VOCS
#relative amounts of visits from pollinator groups
distance.matrix.pollgroup.relative <- vegdist(decostand(net.grp, method="tot"))
mantel(distance.matrix, distance.matrix.pollgroup.relative)
Mantel statistic based on Pearson's product-moment correlation
Call:
mantel(xdis = distance.matrix, ydis = distance.matrix.pollgroup.relative)
Mantel statistic r: -0.1527
Significance: 0.983
Upper quantiles of permutations (null model):
90% 95% 97.5% 99%
0.106 0.142 0.167 0.206
Permutation: free
Number of permutations: 999
#similarity of the relative amounts of scent of each plant
distance.matrix.relative <- vegdist(decostand(scent.net.cut, method="tot"))
#is similarity of the relative amounts of scent of each plant correlated with the similarity in relative visits from each pollinator group
mantel(distance.matrix.relative, distance.matrix.pollgroup.relative)
Mantel statistic based on Pearson's product-moment correlation
Call:
mantel(xdis = distance.matrix.relative, ydis = distance.matrix.pollgroup.relative)
Mantel statistic r: -0.1544
Significance: 0.982
Upper quantiles of permutations (null model):
90% 95% 97.5% 99%
0.109 0.149 0.177 0.210
Permutation: free
Number of permutations: 999
#No correlation between scent similarity and similarity of the proportion of visits by pollinators or pollinator group (p> 0.1)
#not grouped by pollinator class but still relative
distance.matrix.pp.relative <- vegdist(decostand(net.cut, method="tot"))
mantel(distance.matrix.relative, distance.matrix.pp.relative)
Mantel statistic based on Pearson's product-moment correlation
Call:
mantel(xdis = distance.matrix.relative, ydis = distance.matrix.pp.relative)
Mantel statistic r: 0.05562
Significance: 0.181
Upper quantiles of permutations (null model):
90% 95% 97.5% 99%
0.0767 0.1018 0.1208 0.1549
Permutation: free
Number of permutations: 999
#interaction web plot
#plant-voc bipartite
plotweb(sqrt(scent.net), text.rot=90, method="cca", col.high=voc.class.colors[compound.class$major_class], bor.col.high=voc.class.colors[compound.class$major_class], col.low=pcols[plants.majorpoll], bor.col.low=pcols[plants.majorpoll], col.interaction=pal[5], bor.col.interaction=pal[5])

scent.net.cut <- scent.net[rownames(scent.net)%in% paulsplants, ]
scent.mds <- metaMDS(sqrt(scent.net.cut), autotransform = F, trace=F)
scent.op <- ordiplot(scent.mds, type="n")
points(scent.op, what="species", pch=3, col="grey50")
points(scent.op, what="sites", cex=1.5, pch=19, col=pcols[plants.majorpoll])
text(scent.op, what="sites", cex=0.8, col=pcols[plants.majorpoll])

#are scents of plants with different major pollinators different? no...
#do phylogeny, plants in same genus will cluster together need to account for phylogeny
scent.cap <- capscale(sqrt(scent.net.cut)~plants.majorpoll.cut)
anova(scent.cap)
Permutation test for capscale under reduced model
Permutation: free
Number of permutations: 999
Model: capscale(formula = sqrt(scent.net.cut) ~ plants.majorpoll.cut)
Df Variance F Pr(>F)
Model 3 808879 0.5422 0.812
Residual 39 19393005
plot(scent.cap, type="text")

library(indicspecies)
plants.majorpoll.cut <- plants.majorpoll[rownames(scent.net.cut)]
mp <- multipatt(as.data.frame(scent.net.cut), plants.majorpoll.cut, control=how(nperm=999))
summary(mp)
Multilevel pattern analysis
---------------------------
Association function: IndVal.g
Significance level (alpha): 0.05
Total number of species: 354
Selected number of species: 9
Number of species associated to 1 group: 4
Number of species associated to 2 groups: 3
Number of species associated to 3 groups: 2
List of species associated to each combination:
Group butterfly #sps. 1
stat p.value
Carbamic acid, N-phenyl-, 1,5-dimethyl-1-vinyl-4-hexenyl ester 0.942 0.035 *
Group hummingbird #sps. 3
stat p.value
Propanoic acid, 2-methyl-, 3-methylbutyl ester 1.000 0.043 *
Petasitene 0.996 0.043 *
2H-Pyran-2-one, tetrahydro-3,6-dimethyl- 0.962 0.037 *
Group butterfly+hummingbird #sps. 2
stat p.value
Bicyclo[7.2.0]undec-4-ene, 4,11,11-trimethyl-8-methylene- 0.967 0.001 ***
Methoxyacetic acid, tetradecyl ester 0.941 0.019 *
Group fly+hummingbird #sps. 1
stat p.value
2-Butenal, 3-methyl- 0.921 0.045 *
Group butterfly+fly+hummingbird #sps. 2
stat p.value
(E)-4-Oxohex-2-enal 0.941 0.008 **
1-Propanol, 2-methyl- 0.890 0.030 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
###Do generalist plants have higher diversity of scent compounds?
plants.numlinks <- rowSums(decostand(net.cut, "pa"))
#polls.classtree <- compute.brlen(polls.classtree$phylo,1)
#library(picante)
#plants.pdlinks <- pd(net.lump, polls.classtree)$PD
plants.numcompounds <- rowSums(decostand(scent.net.cut, "pa"))
plants.sumcompounds <- rowSums(scent.net.cut)
plants.shannoncompounds <- diversity(scent.net.cut, index="shannon")
plants.shannonlinks <- diversity(net.grp, index="shannon")
scent.diversity <- tibble(numlinks=plants.numlinks, numcompounds=plants.numcompounds, sumcompounds=plants.sumcompounds, shannoncompounds=plants.shannoncompounds, shannonlinks=plants.shannonlinks)
ggplot(scent.diversity, aes(y=numcompounds, x=numlinks))+
geom_point()+ geom_smooth()

anova(lm(numlinks~numcompounds, data=scent.diversity))
Analysis of Variance Table
Response: numlinks
Df Sum Sq Mean Sq F value Pr(>F)
numcompounds 1 792.3 792.30 3.8827 0.05556 .
Residuals 41 8366.4 204.06
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#tells us if plants with more diverse scent get more pollinator species
anova(lm(plants.numlinks~shannoncompounds, data=scent.diversity))
Analysis of Variance Table
Response: plants.numlinks
Df Sum Sq Mean Sq F value Pr(>F)
shannoncompounds 1 35.4 35.438 0.1593 0.6919
Residuals 41 9123.3 222.520
#plot(plants.numcompounds~plants.pdlinks)
#plot(plants.shannoncompounds~plants.pdlinks)
shannon.compoundxlinks.plot <- ggplot(scent.diversity, aes(y=shannonlinks, x=shannoncompounds))+
geom_point()+ geom_smooth()
ggsave(shannon.compoundxlinks.plot, file='shannon.compoundxlinks.plot.jpg', width=130, height=105, units= 'mm', dpi=300)
compoundxlinks.plot <- ggplot(scent.diversity, aes(y=plants.numlinks, x=shannoncompounds))+
geom_point()+ geom_smooth()
ggsave(compoundxlinks.plot, file='compoundxlinks.plot.jpg', width=130, height=105, units= 'mm', dpi=300)
shannon.compound.plot <- ggplot(scent.diversity, aes(y=shannoncompounds, x=plants.majorpoll.cut))+
geom_point()+ geom_smooth()
ggsave(shannon.compound.plot, file='shannon.compound.plot.jpg', width=110, height=105, units= 'mm', dpi=300)
numcompounds.plot <-ggplot(scent.diversity, aes(y=numcompounds, x=plants.majorpoll.cut))+
geom_point()+ geom_smooth()
ggsave(numcompounds.plot, file='numcompounds.plot.jpg', width=110, height=105, units= 'mm', dpi=300)
sumcompounds.plot<- ggplot(scent.diversity, aes(y=sumcompounds, x=plants.majorpoll.cut))+
geom_point()+ geom_smooth()
ggsave(sumcompounds.plot, file='sumcompounds.plot.jpg', width=110, height=105, units= 'mm', dpi=300)
#interaction web plot
#voc-pollinator bipartite
#make dataframe with vocs and pollinators
##change scent.net to long format
plantcompoundemission <- scent.net.cut %>% as.data.frame() %>% rownames_to_column("plant") %>% #average peak area
pivot_longer(-plant, names_to = "compound", values_to = "emissions") %>% filter(emissions!=0)
#write_csv(plantcompoundemission, "plantcompoundemissions.csv")
#total number of visits to plant by each pollinator species
plantpollinatorvisits <- net.cut %>% as.data.frame() %>% rownames_to_column("plant") %>% pivot_longer(-plant, names_to = "pollinator", values_to = "visits")
#total number of visits to each compound. Pollinators rows, compounds cols, values = number of visits
scentpollnet <- left_join(plantcompoundemission, plantpollinatorvisits,relationship = "many-to-many") %>% group_by(pollinator, compound) %>% summarise(visits=sum(visits)) %>% pivot_wider(names_from = "compound", values_from = "visits") %>% column_to_rownames("pollinator") %>% as.matrix()
#relative amount of each compound
plantcompoundrelative <- scent.net.cut %>% as.data.frame() %>% decostand(method="tot") %>% rownames_to_column("plant") %>% #average peak area
pivot_longer(-plant, names_to = "compound", values_to = "emissions") %>% filter(emissions!=0)
#multiply number of visits by relative peak area
# weights visit to compounds. Proportions visits based on how many VOCs there are.
scentpollnet.area <- left_join(plantcompoundrelative, plantpollinatorvisits,relationship = "many-to-many") %>% group_by(pollinator, compound) %>% mutate(visits.area=visits*emissions) %>% summarise(visits.area=sum(visits.area)) %>% pivot_wider(names_from = "compound", values_from = "visits.area") %>% column_to_rownames("pollinator") %>% as.matrix()
scentpollnet.class <- tibble(pollinator=rownames(scentpollnet)) %>% left_join(polls.class.dataframe) %>% pull(pollinator_group)
#plot bipartite network
#plotweb(sqrt(scentpollnet), text.rot=90, method="cca", col.high=pal[2], bor.col.high=pal[2], #col.low=pcols[scentpollnet.class], bor.col.low=pcols[scentpollnet.class], col.interaction=pal[5], #bor.col.interaction=pal[5])
#plotweb for scentpollnet.area
plotweb(sqrt(scentpollnet.area), text.rot=90, method="cca", col.high=pal[2], bor.col.high=pal[2], col.low=pcols[scentpollnet.class], bor.col.low=pcols[scentpollnet.class], col.interaction=pal[5], bor.col.interaction=pal[5])

#value = number of
heatmaply(scentpollnet, scale="none", showticklabels = c(F,T), margins= c(0,NA,0,0),
scale_fill_gradient_fun = ggplot2::scale_fill_gradient(low = "white", high = "#000050", trans="log1p"),
distfun=vegdist,
hclustfun=function(x) hclust(x, method="mcquitty"),
#row_side_colors = data.frame(P=plants.majorpoll[rowSums(net)>0]),
#row_side_palette=function(n) pcols,
row_side_colors=data.frame(PollinatorGroup=scentpollnet.class),
row_side_palette=function(n) pcols.all)
voc.pollgrp.cap <- capscale(sqrt(scentpollnet.area)~scentpollnet.class)
anova(voc.pollgrp.cap)
Permutation test for capscale under reduced model
Permutation: free
Number of permutations: 999
Model: capscale(formula = sqrt(scentpollnet.area) ~ scentpollnet.class)
Df Variance F Pr(>F)
Model 4 15.904 1.67 0.111
Residual 88 209.505
plot(voc.pollgrp.cap, type="text")

#make vector with compound class and all scent compounds
vol.class <-deframe(select(compound.class, compound_long, major_class) )
#turns scentpollnet matrix to data frame then pivots longer to have 2 cols called chemical and visits
chemical.taxa <- left_join(plantcompoundemission, plantpollinatorvisits,relationship = "many-to-many") %>% left_join(polls.class.dataframe) %>% left_join(compound.class, by=c("compound"= "compound_long")) %>% group_by(pollinator_group, major_class) %>% summarise(visits=sum(visits)) %>% pivot_wider(names_from = major_class, values_from = visits) %>% column_to_rownames("pollinator_group") %>% as.matrix()
#chemical class & pollinator group network
# how do I incorporate chemical class instead of compounds?
#make data frame with chemical class and pollinator family - use scent.net & chemical.class
#and use scentpollnet but scentpollnet compounds in different order than scent.net
plotweb(sqrt(chemical.taxa), text.rot=90, method="normal", col.high=pal[3], bor.col.high=pal[3], col.low=pcols.all, bor.col.low=pcols.all, col.interaction=pal[5], bor.col.interaction=pal[5])

library(bipartite)
library(vcd)
#side quest
mosaic(chemical.taxa)

mosaic(net)

chisq.test(chemical.taxa) #p=1 so this means that no specialist compound classes. and vise versa. nobody is specialized
Pearson's Chi-squared test
data: chemical.taxa
X-squared = 8488.8, df = 24, p-value < 2.2e-16
#value = number of
heatmaply(scentpollnet.area, scale="none", showticklabels = c(F,T), margins= c(0,NA,0,0),
scale_fill_gradient_fun = ggplot2::scale_fill_gradient(low = "white", high = "#000050", trans="log1p"),
distfun=vegdist,
hclustfun=function(x) hclust(x, method="mcquitty"),
#row_side_colors = data.frame(P=plants.majorpoll[rowSums(net)>0]),
#row_side_palette=function(n) pcols,
row_side_colors=data.frame(PollinatorGroup=scentpollnet.class),
row_side_palette=function(n) pcols.all)
#species accumulation curves
#figure out what percent Pauls data catpured the whole pollinator community
#pull out 1 random 1 hr sampling period and say how many pollinators observered, then do this again and again randomly - variation in samples = standard error y axis= number of species, x= number of pollinator observations
#class= bee, fly, butterfly, hummingbird
#xaxis is number of plant species observed (if you did the whole summer but only looked at 1 species and then next summer looked at 2 species, etc.)
#
#simplify net.long to only have numbers in it so that accumulation curve will work
#rows are observation "units" and columns are pollinators
net.wide <- net.long %>% select(pollinator, interactions) %>% mutate(index=row_number()) %>%
pivot_wider(names_from=pollinator, values_from = interactions) %>% mutate(across(everything(),~replace_na(.x, 0))) %>% select(-index) %>% as.matrix()
for(class in unique(polls.class)) { #could replace class with family
class_subset <- net[,polls.class==class] #making accumulation curve for each pollinator class
class_sa <- vegan::specaccum(class_subset, permutations = 999)
plot.new()
plot(class_sa, col = as.integer(factor(class))+1, add = class != "bee", lwd=2,
xlab="Pollinator Observations", ylab="Pollinator Species", main="Pollinator Species accumulation curve", sub=class)
}




legend("bottomright", unique(polls.class),
col = 2:5, pch = 19, title = "Class")

sa<- vegan::specaccum(net.wide,permutations=100, method="exact")
plot(sa)

#Pauls total pollinator observations 8hrs/wk * 42wks = 336 hours total
visitationrate <- nrow(net.long)/336 #15.33 unique p-p observations per hr (observation can include multiple visits)
#observation=within a 15min period how many p-p interactions did we see, total number of links observed
2000/visitationrate # 130 hrs per treatment ( 2 treatments = 260 hrs per year)
[1] 130.4854
#then based on my observation hours, what percent of the pollinator community could be captured
library(coin)
library(rstatix)
ND_poll<-ND(net, normalised = T)$higher
#creates dataframe with normalized degree for each pollinator and pollinator class (same as data.frame)
ND.poll.class <- tibble(ND=ND_poll, polls.class=polls.class)
ND.poll.class$ND <- as.numeric(ND.poll.class$ND)
ND.poll.class$polls.class <- as.factor(ND.poll.class$polls.class)
kruskal_test(ND~polls.class, data=ND.poll.class)
Asymptotic Kruskal-Wallis Test
data: ND by
polls.class (bee, butterfly, fly, hummingbird, wasp)
chi-squared = 4.0146, df = 4, p-value = 0.404
ggplot(ND.poll.class, aes(x= polls.class, y=ND)) +geom_boxplot()

library(rstatix)
#ND for pollinators in scent network
ND_poll_scent<-ND(scentpollnet, normalised = T)$lower
#ND of VOC chemical class
ND_scent_poll<-ND(scentpollnet, normalised = T)$higher
#creates dataframe with normalized degree for each pollinator and pollinator class (same as data.frame)
ND.poll <- tibble(ND=ND_poll_scent, polls.class=polls.class)
#kruskal_test doesn't like characters
ND.poll$ND <- as.numeric(ND.poll$ND)
ND.poll$polls.class <- as.factor(ND.poll$polls.class)
kruskal_test(ND~polls.class, data=ND.poll)
Asymptotic Kruskal-Wallis Test
data: ND by
polls.class (bee, butterfly, fly, hummingbird, wasp)
chi-squared = 2.6914, df = 4, p-value = 0.6107
ggplot(ND.poll, aes(x= polls.class, y=ND)) +geom_boxplot()

#creates dataframe with normalized degree for each VOC and chemical class (same as data.frame)
ND.scent <- tibble(ND=ND_scent_poll, vols.class=compound.class$major_class)
kruskal_test(ND~vols.class, data=ND.scent)
Asymptotic Kruskal-Wallis Test
data: ND by
vols.class (aliphatic, monoterpene, sesquiterpene, nitrogenous, sulfur, benzenoid, other)
chi-squared = 7.8313, df = 6, p-value = 0.2507
ggplot(ND.scent, aes(x= vols.class, y=ND)) +geom_boxplot()

#created a function similar to normalized degree(ND) but ND looks at presence/absence in the network where normalized visits (NV) looks at the normalized weighted visits/interactions. This accounts for variation in the amount each plant or compound is visited and give those with higher interaction frequency greater weight or contribution to the network
normalized.visits <- function (web, normalised = TRUE)
{
websum <- sum(web) #takes sum of all interactions in the web
dlower <- rowSums(web)
dhigher <- colSums(web)
Nlow <- Nhigh <- 2
if (normalised) {
Nlow <- websum
Nhigh <- websum
}
low <- dlower/Nlow #weighting the interactions by dividing rowsums by the sum of all links in the network
names(low) <- rownames(web) #dividing colsums by the sum of all links in the network
high <- dhigher/Nhigh
names(high) <- colnames(web)
list(lower = low, higher = high)
}
#Which pollinators have the most visits. of all the visits recorded, what percentage of them were by that pollinator (for a given plant what percentage of visits does it get out of the whole meadow, for a pollinator what percentage of visits does it give out of all visits)
#normalized visits for pollinators in scent net
NV_poll<-normalized.visits(net, normalised = T)$higher
#creates dataframe with normalized degree for each pollinator and pollinator class (same as data.frame)
NV.poll.class <- tibble(NV=NV_poll, polls.class=polls.class)
#kruskal_test won't work as.character, need to change to numeric and factor
NV.poll.class$NV <- as.numeric(NV.poll.class$NV)
NV.poll.class$polls.class <- as.factor(NV.poll.class$polls.class)
kruskal_test(NV~polls.class, data=NV.poll.class)
Asymptotic Kruskal-Wallis Test
data: NV by
polls.class (bee, butterfly, fly, hummingbird, wasp)
chi-squared = 11.58, df = 4, p-value = 0.02076
ggplot(NV.poll.class, aes(x= polls.class, y=NV)) +geom_boxplot()
-1.png)
#only a few species making most of the visits
#1 bee responsible for over 40% of visits
#normalized visits for pollinators in scent net
#which pollinator class visits the most amount of compounds weighted by the number of visits
NV_poll_scentnet<-normalized.visits(scentpollnet, normalised = T)$lower
#creates dataframe with normalized degree for each pollinator and pollinator class (same as data.frame)
NV.poll.class.scentnet <- tibble(NV=NV_poll_scentnet, polls.class=polls.class)
#kruskal_test won't work as.character, need to change to numeric and factor
NV.poll.class.scentnet$NV <- as.numeric(NV.poll.class.scentnet$NV)
NV.poll.class.scentnet$polls.class <- as.factor(NV.poll.class.scentnet$polls.class)
kruskal_test(NV~polls.class, data= NV.poll.class.scentnet)
Asymptotic Kruskal-Wallis Test
data: NV by
polls.class (bee, butterfly, fly, hummingbird, wasp)
chi-squared = 2.9103, df = 4, p-value = 0.573
ggplot(NV.poll.class.scentnet, aes(x= polls.class, y=NV)) +geom_boxplot()

#super generalist bees
#normalized visits for VOCs in scent net
#which compound class are most visited
NV_voc<-normalized.visits(scentpollnet, normalised = T)$higher
#creates dataframe with normalized degree for each pollinator and pollinator class (same as data.frame)
NV.vols.class <- tibble(NV=NV_voc, vols.class=compound.class$major_class)
k_t.vols.class <- kruskal_test(NV~vols.class, data=NV.vols.class)
library(coin)
library(PMCMRplus) #for posthoc
library(PMCMR)
#posthoc_vols.class <- posthoc.kruskal.nemenyi.test(NV.vols.class$vols.class, NV.vols.class$NV, dist = "tukey")
ggplot(NV.vols.class, aes(x= vols.class, y=NV)) +geom_boxplot()

#tells you what vocs get the most visits
#no difference in visitation rate to compound classes
#Modularity
library(vcd)
modules.net <- computeModules(sqrt(net))
plotModuleWeb(modules.net, labsize=.5)

module.net.list <- listModuleInformation(modules.net)
flatten(module.net.list[[2]])
[[1]]
[1] "androsace_septentrionalis" "claytonia_lanceolata"
[3] "draba_aurea" "erythronium_grandiflorum"
[5] "hydrophyllum_capitatum" "lomatium_dissectum"
[7] "mertensia_fusiformis" "pseudocymopterus_montanus"
[9] "ranunculus_inamoenus" "taraxacum_officinale"
[11] "viola_praemorsa" "amelanchier_alnifolia"
[13] "fragaria_virginiana" "hydrophyllum_fendleri"
[15] "linum_lewisii"
[[2]]
[1] "andrena_cyanophila" "lasioglossum_spp" "sphecodes_sp_01"
[4] "syrphidae_spp" "osmia_iridis" "scathophagidae_sp_01"
[7] "halictus_confusus" "halictus_virgatellus" "halictidae_green_spp_1"
[10] "eupeodes_volucris" "andrena_sp_02" "muscidae_sp_02"
[13] "syrphidae_02" "rhamphomyia_sp_01" "cartosyrphus_tarda"
[16] "callophrys_affinis" "osmia_lignaria" "cartosyrphus_sp_02"
[19] "tephritidae_01" "polygonia_zephyrus"
[[3]]
[1] "delphinium_nuttallianum" "ipomopsis_aggregata"
[3] "lathyrus_lanszwertii" "solidago_multiradiata"
[5] "gentiana_parryi" "vicia_americana"
[[4]]
[1] "bombus_appositus" "selasphorus_platycercus"
[3] "coelioxus_funeraria" "anthidium_emarginatum"
[5] "bombus_nevadensis" "papilio_gothica"
[7] "sphex_spp_1"
[[5]]
[1] "agoseris_aurantiaca" "campanula_rotundifolia"
[3] "erigeron_speciosus" "helianthella_quinquenervis"
[5] "heliomeris_multiflora" "heterotheca_villosa"
[7] "senecio_serra" "castilleja_sulphurea"
[9] "phacelia_heterophylla"
[[6]]
[1] "nymphalidae_spp" "bombus_bifarius"
[3] "megachile_relativa" "tachinidae_sp_01"
[5] "bombus_sylvicola" "lycaenidae_spp"
[7] "anthophora_terminalis" "bombus_flavifrons"
[9] "andrena_costillensis" "bombus_californicus"
[11] "bombus_occidentalis" "megachile_pugnata"
[13] "arctophila_flagrans" "osmia_grindeliae"
[15] "pieridae_spp" "speyeria_mormonia"
[17] "chrysotoxum_ventricosum" "megachile_melanophea"
[19] "psithyrus_insularis" "lycaena_rubidus_sirius"
[21] "colias_alexandra" "bombus_frigidus"
[23] "mesebrina_latreillei" "syrphidae_spp_2"
[25] "syrphidae_spp_4" "bombus_balteatus"
[27] "syrphidae_spp_1" "phormia_spp_1"
[29] "megachile_frigida" "lycaena_cupreus"
[31] "bombus_mixtus" "osmia_tristela"
[33] "melissodes_confusa"
[[7]]
[1] "agoseris_glauca" "arenaria_congesta" "eriogonum_subalpinum"
[4] "eriogonum_umbellatum" "mahonia_repens" "mertensia_ciliata"
[7] "potentilla_hippiana" "rosa_woodsii" "sedum_lanceolatum"
[10] "achillea_millefolium" "potentilla_gracilis"
[[8]]
[1] "muscidae_spp_01" "colletes_kincaidii"
[3] "halictus_rubicundus" "andrena_transnigra"
[5] "halictidae_spp_01" "panurginus_ineptus"
[7] "thricops_septentrionalis" "osmia_bucephala"
[9] "chrysis_coerulans" "villa_eumenes"
[11] "glaucopsyche_lygdamus" "pieris_mcdunnoughi"
[13] "andrena_vicinoides" "colletes_nigrifrons"
[15] "melanstoma_caerulescens"
[[9]]
[1] "boechera_stricta" "erigeron_flagellaris" "galium_boreale"
[4] "senecio_integerrimus"
[[10]]
[1] "ochlodes_sylvanoides" "systoechus_sp" "osmia_proxima"
[4] "pieris_rapa" "hoplitus_fulgida" "osmia_subaustralis"
[7] "osmia_montana" "osmia_coloradensis" "sphaerophoria_robusta"
[10] "melanostoma_kelloggi" "syrphidae_spp_3" "eupeodes_lapponicus"
[13] "euphydras_anicia" "eristalis_latifrons" "osmia_sp_sgo"
[16] "cryptopogon_sp_01" "hoplitus_robusta" "coenonympha_ochracea"
pollinator.modules <- tibble(module=1:5) %>% mutate(pollinator=map(module, ~module.net.list[[2]][[.x]][[2]])) %>% unnest(pollinator) %>% left_join(polls.class.dataframe) %>% count(module, pollinator_group) %>% pivot_wider(names_from = pollinator_group, values_from = n, values_fill = 0) %>% column_to_rownames("module") %>% as.matrix()
chisqr.polls <- chisq.test(pollinator.modules)
residules.chisqr.polls <- residuals(chisqr.polls)
mosaic(t(pollinator.modules), shade = T, gp=shading_hcl)

#test to see p-p network modules
#mosaic(net, shade = T, gp=shading_hcl)
modules.scentnet <- computeModules(sqrt(t(scentpollnet)))
#plotModuleWeb(modules.scentnet)
module.scentnet.list <- listModuleInformation(modules.scentnet)
pollinator.modules <- tibble(module=1:4) %>% mutate(pollinator=map(module, ~module.scentnet.list[[2]][[.x]][[2]])) %>% unnest(pollinator) %>% left_join(polls.class.dataframe) %>% count(module, pollinator_group) %>% pivot_wider(names_from = pollinator_group, values_from = n, values_fill = 0) %>% column_to_rownames("module") %>% as.matrix()
chisqr.polls <- chisq.test(pollinator.modules)
residules.chisqr.polls <- residuals(chisqr.polls)
mosaic(t(pollinator.modules), shade = T, gp=shading_hcl, labeling = vcd::labeling_border(rot_labels = c(0, 0)))

#chemical class modularity
chemical.modules <- tibble(module=1:4) %>% mutate(compound_long=map(module, ~module.scentnet.list[[2]][[.x]][[1]])) %>% unnest(compound_long) %>% left_join(compound.class) %>% count(module, major_class) %>% pivot_wider(names_from = major_class, values_from = n, values_fill = 0) %>% column_to_rownames("module") %>% as.matrix()
mosaic(t(chemical.modules), shade = T, gp=shading_hcl, labeling = vcd::labeling_border(rot_labels = c(0, 0)))

#Stacked bar chart with flowers visited by each pollinator #Looks at the top 20 plants and top 75 pollinators (interactions square rooted to adjust scale on the figure) #x= interactions, y=pollinators, color=plant
#Stacked bar chart with what pollinators visited each plant species #45 plant species, 18 pollinators #x= interactions, y=plant, color=pollinator
library(ggplot2)
library(tidyverse)
library(RColorBrewer)
mypalette= c(brewer.pal(8, "Set2"), brewer.pal(12, "Set3"))
#looks at top 20 plants and top 75 pollinators (interactions sqrt to adjust scale on figure)
#pollinator perspective
ggplot(net.long %>% mutate(plant=fct_lump_n(plant, 18), pollinator=fct_lump_n(pollinator, 75)), aes(x=sqrt(interactions), y=pollinator, fill=plant)) + geom_col() + scale_fill_manual(values=mypalette)

#plant perspective
ggplot(net.long %>% mutate(plant=fct_lump_n(plant, 45), pollinator=fct_lump_n(pollinator, 18)), aes(x=sqrt(interactions), y=plant, fill=pollinator)) + geom_col() + scale_fill_manual(values=mypalette)

#Uses Paul’s phenology data to create stackplot
#Panels are years, x= week number, y= flower count, color= plant species ## Plot 2 - flower phenology heatmap ## Plot 3 - Pollinator Phenology #Panels are years, x= week number, y= interactions, color= pollinator ## Plot 4 - flower phenology heatmap ## Plot 5 - Scent Phenology #weighted average of flower scents and mean is weighted by how many flowers there are that week #shows relative amount of that compound out of the total scent in the meadow that week #each color is a compound, value= total emissions ## Plot 6 - heatmap scent phenology
phenology <- read_sheet("15Jnx8OWMzpakC0YvJWmZMmaRSnA0X-7HBjymDsvUfjY")
#stack plots and heatmaps
#Flowers#
ggplot(phenology, aes(x=week_num, y=flower_count, fill= plant )) + geom_col() + scale_fill_manual(values=sample(rainbow(46))) +facet_wrap(vars(year))

#heatmap
ggplot(phenology, aes(x=week_num, y=plant, fill= sqrt(flower_count ))) + geom_tile() + scale_fill_viridis_c() + facet_wrap(vars(year))

#pollinators#
ggplot(net.long, aes(x=week_num, y=interactions, fill= pollinator )) + geom_col() + scale_fill_manual(values=sample(rainbow(93))) +facet_wrap(vars(year))

#heatmap
ggplot(net.long, aes(x=week_num, y=pollinator, fill= sqrt(interactions ))) + geom_tile() + scale_fill_viridis_c() + facet_wrap(vars(year))

#scent#
#number of flowers per week and multiple by mean scent of that species
#how many flowers are open each week
flowers.per.week <- phenology %>% group_by(year, week_num, plant) %>%
summarise(flower_count=sum(flower_count))
#ideally want for each week, 1 average weighted mean of scent compound
#but we only have a few scent samples so just taking the mean of scent compounds for each plant
#flowers per week percentage
flower.percentage <- flowers.per.week %>% group_by(year, week_num) %>%
mutate(flower_count=flower_count/sum(flower_count)) %>%
pivot_wider(names_from = c(year,week_num), values_from = flower_count, values_fill = 0) %>%
filter(plant%in%rownames(scent.net)) %>%
arrange(match(plant, rownames(scent.net)))%>% #matches scent.net plant order
column_to_rownames("plant") %>% as.matrix()
#change names in pauls data
scentxphenology <- t(flower.percentage) %*% scent.net[rownames(flower.percentage),]
#values in matrix = weighted average of flower scents and mean is weighted by how many flowers there are that week
#replaces Nans with 0
scentxphenology[is.nan(scentxphenology)] <- 0
#ggplot likes long format
scent.phenology <- scentxphenology %>% as.data.frame() %>% rownames_to_column("yearweek") %>% separate(yearweek, into = c("year", "week_num")) %>% mutate(year=as.integer(year), week_num= as.integer(week_num)) %>% pivot_longer(all_of(colnames(scent.net)))
ggplot(scent.phenology, aes(x=week_num, y=value, fill= name )) + geom_col() + scale_fill_manual(values=sample(rainbow(355)), guide="none") +facet_wrap(vars(year))

#heatmap - scent x phenology
ggplot(scent.phenology, aes(x=week_num, y=name, fill= sqrt(value ))) + geom_tile() + scale_fill_viridis_c() + facet_wrap(vars(year))

#number of compounds each plant species has
#total emissions for each species
#diversity of scents (shannon diversity)
library(tidyverse)
scent.net.long <- scent.net %>% as.data.frame() %>% rownames_to_column("plant") %>% #average peak area
pivot_longer(-plant, names_to = "compound", values_to = "emissions") %>% filter(emissions!=0)
compound.class.percentage <- scent.net.long %>% left_join(compound.class, by=c(compound="compound_long")) %>% group_by(plant, major_class) %>% summarise(emissions=sum(emissions)) %>% #gives total emissions for each compound class
mutate(emissions=emissions/sum(emissions)) %>% #relative amount of each compound class in species scent
pivot_wider(names_from = major_class, values_from = emissions, names_prefix = "percent_", values_fill = 0) #makes column for each compound class and gives percentage of that compound class out of the total emissions for that species
#rare vs common compounds (mean rarity)
##for each compound, how many times does it occur in the dataset?
compound.rarity <- scent.net.long %>% count(compound) %>% #number of species the compound occurs in
mutate(rarity=n/nrow(scent.net)) #the percentage of species the compound occurs in
weighted.scent.rarity <- scent.net.long %>% left_join(compound.rarity) %>% group_by(plant) %>%
mutate(emissions=emissions/sum(emissions)) %>% #makes emissions -> relative emissions
mutate(weighted_scent_rarity = emissions*rarity) %>% #multiplies emissions x rarity to get weighted average
summarise(weighted_scent_rarity= sum(weighted_scent_rarity)) #then add the different parts of the average
#rarity score of 1 means this compound is found in all the plant species, 0= no other species has these compounds
#mertensia has 0.70 rarity - the weighted average of all the compounds rarity; 70% of compounds (or total emissions) are shared, 30% of compounds or total emissions other species don't have
scent.summary <- tibble(plant=rownames(scent.net), total_emissions= rowSums(scent.net), number_compounds=rowSums(scent.net>0),scent_shannon_diversity= diversity(scent.net, index="shannon")) %>%
left_join(weighted.scent.rarity) %>%
left_join(compound.class.percentage)
#pollinators that visit plants
#Percentage of visits by each pollinator family to each plant species
pollinator.family.percentage <- net.long %>% group_by(plant, pollinator_family) %>% summarise(interactions=sum(interactions)) %>% #gives sum of interactions for each pollinator family
mutate(interactions=interactions/sum(interactions)) %>% #relative amount of each pollinator family visiting each flower species
pivot_wider(names_from = pollinator_family, values_from = interactions, names_prefix = "percent_", values_fill = 0)
#Percentage of visits by each pollinator group to each plant species
pollinator.group.percentage <- net.long %>% group_by(plant, pollinator_group) %>% summarise(interactions=sum(interactions)) %>% #gives sum of interactions for each pollinator group
mutate(interactions=interactions/sum(interactions)) %>% #relative amount of each pollinator group visiting each flower species
pivot_wider(names_from = pollinator_group, values_from = interactions, names_prefix = "percent_", values_fill = 0)
#rarity - how generalized the pollinator is
pollinator.rarity <- net.long %>% count(plant, pollinator) %>% #adds how many rows of interactions between plant & poll species
count(pollinator, name="number_plants_visited") %>% #number of plants the pollinator goes to
mutate(rarity=number_plants_visited/nrow(net)) #the percentage of plant species the pollinator visits in the meadow
#gives for each plant species, which pollinators visited and their relative contributions to total number of visits
net.long.relative.interactions <- net.long %>% group_by(plant, pollinator) %>%
summarise(interactions=sum(interactions)) %>%
left_join(pollinator.rarity) %>%
group_by(plant) %>%
mutate(interactions=interactions/sum(interactions))
#percent of interactions that come from each pollinator species
weighted.pollinator.rarity <- net.long.relative.interactions %>% mutate(weighted_poll_rarity = interactions*rarity) %>% #multiplies interactions x rarity to get weighted average
summarise(weighted_poll_rarity = sum(weighted_poll_rarity))
flower.visitor.summary <- tibble(plant=rownames(net), total_visits= rowSums(net), number_pollspecies_visited=rowSums(net>0),pollinator_shannon_diversity= diversity(net, index="shannon")) %>%
left_join(weighted.pollinator.rarity) %>%
left_join(pollinator.family.percentage) %>% left_join(pollinator.group.percentage)
megatable <- scent.summary %>% full_join(flower.visitor.summary) #full join bc pauls data and scent data plants don't match
library(pheatmap)
library(smatr)
correlation.table <- megatable %>% column_to_rownames("plant") %>% cor(use="na.or.complete")
pheatmap(correlation.table, scale="none", clustering_method = "mcquitty", color=colorRampPalette(c("blue","white", "red"))(200),
breaks = seq(-1,1, by=0.01)) #which ones it labels on the scale

#scatter plot with just weighted scent and weighted pollinator visits
weighted.scent.poll <- ggplot(megatable, aes(x=weighted_poll_rarity, y=weighted_scent_rarity, color=number_pollspecies_visited)) + geom_point()+geom_smooth() + geom_text(aes(label=plant))
print(weighted.scent.poll)

#the most common volatiles get the most visits from the most generalized pollinators
sma.weighted.scent.poll <- sma(weighted_poll_rarity~weighted_scent_rarity, data=megatable )
summary(sma.weighted.scent.poll)
Call: sma(formula = weighted_poll_rarity ~ weighted_scent_rarity, data = megatable)
Fit using Standardized Major Axis
------------------------------------------------------------
Coefficients:
elevation slope
estimate -0.3602556 1.352232
lower limit -0.5885207 1.001742
upper limit -0.1319906 1.825352
H0 : variables uncorrelated
R-squared : 0.06769348
P-value : 0.091989
#percent_hesperiidae x nitrogenous vocs
ggplot(megatable, aes(x=percent_nitrogenous, y=percent_lycaenidae, color=number_pollspecies_visited)) + geom_point()+geom_smooth() + geom_text(aes(label=plant))

##pollinator trait summary